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Abstract 

{Chung, 2009 #1}The transforming growth factor-p (TGF-P) superfamily of cytokines plays a fundamental role in a wide 
variety of cellular processes, including growth, differentiation, apoptosis, and tissue homeostasis. Its relevance is 
emphasized by the mutations of its core components that are associated with diverse human diseases, such as cancer and 
cardiovascular pathologies. A prominent regulator of the pathway is Smad7, which attenuates the signal and controls its 
duration in a cell-type-dependent manner through a negative feedback loop. Here, we characterize all the potential Smad7- 
mediated negative feedback network motifs and investigate their effects on the signaling dynamics upon stimulation with 
TGF-p and bone morphogenetic protein (BMP) ligands. The results show that the specific negative feedback 
implementation is a key determinant of both the response of the system to single and multiple ligands of the TGF-p 
superfamily and its robustness and sensitivity to parameter perturbations. 
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Introduction 

The transforming growth factor- P (TGF-P) superfamily, which 
comprises 33 different ligands in mammalian cells, plays a 
fundamental role in development and maintenance of tissue 
homeostasis [1]. These ligands regulate key cellular processes, such 
as proliferation, motility, differentiation, and apoptosis [2]. 
Dysregulation of the TGF-P signal transduction pathway resulting 
from mutations of its core components has been associated with a 
number of human diseases, including cancer and vascular 
disorders [3,4] . As a result, a significant effort of clinical research 
focuses on developing therapies targeting the TGF-P pathway [1]. 

The multiple cellular effects elicited by the TGF-P superfamily 
ligands are triggered by binding of the ligand to two types of 
receptor serine/threonine protein kinases (type II and type I 
receptors) at the plasma membrane, which then form an active 
ligand-receptor complex. The signal is thenceforth propagated 
through the intracellular Smad proteins into the nucleus where 
activated Smad complexes act as transcription factors, controlling 
the expression of hundreds of genes in a cell-type and context 
dependent way [2]. Specifically, the active ligand-receptor 
complex is internalized into early endosomes, where it recruits 
and phosphorylates one of the receptor-regulated Smad (R-Smad) 
proteins. Phosphorylated R-Smads bind Smad4, the common- 
mediator Smad, forming a heterooligomer that translocates into 
the nucleus and binds to DNA to regulate the expression of its 
target genes. Ligands of the TGF-P superfamily signal through the 
activation of two parallel R-Smad channels. Specifically, bone 
morphogenetic proteins (BMPs) activate the Smad 1/5/8 channel; 
nodal and activin ligands activate the Smad2/3 channel; and 
TGF-P activates both channels [5]. 



Inhibitory Smads (I-Smads), Smad6 and Smad7, negatively 
regulate signaling in this pathway, antagonizing the effects of R- 
Smads [6]. They inhibit the signal through different mechanisms, 
such as sequestering phosphorylated R-Smads, specifically Smadl, 
in an inactive complex as observed for Smad6 [7] or, more 
typically, by competing with R-Smads for receptor binding [8,9] 
and promoting degradation of ligand-receptor complexes through 
Smurf-dependent ubiquitination [10,11]. Importantly, this inhibi- 
tion can occur through a negative feedback loop because TGF-P 
superfamily ligands induce transcription of Smad6 and Smad7 genes 
by the binding of nuclear phosphorylated R-Smad-Smad4 
complexes to the promoters [1]. 

In the last few years, several mathematical models of the Smad- 
dependent TGF-P signal transduction pathway have been 
developed to get insights into its functioning [12-25]. In particular, 
a few of these computational models have incorporated the 
Smad7-mediated negative feedback loop as an explicit component 
of the pathway in order to investigate its mechanistic role in the 
observed behavior [12,15,16,24,25]. In these cases, the effects of 
Smad6, which preferentially blocks BMP signaling, are typically 
combined with those of Smad7, which blocks both TGF-P and 
BMP signaling, in a single effective inhibitory component. We 
have previously demonstrated that differences in the implemen- 
tation of the negative feedback loop capture the distinct signaling 
dynamics of diverse cell lines [25]. In addition to investigating the 
dynamic response of TGF-P signaling, quantitative approaches 
have revealed how the signaling behavior is affected by 
perturbations of its parameters through the use of sensitivity 
analyses [12-15,17,22,25], analytical calculations [18], and other 
types of mathematical analysis [19]. 
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In mammalian cells, the Smad7-dependent negative feedback 
loop has different implementations in different cell lines [25]. For 
example, bovine aortic endothelial cells (BAECs) exhibit an auto- 
regulatory negative feedback loop, where Smad7 is expressed 
through activation of the Smadl channel and inhibits further 
activation of the same R-Smad channel [26]. The mouse myoblast 
cell line C2C12 displays a cross-regulatory negative feedback, 
where Smad7 is expressed through activation of the Smad2 
channel, but inhibits the Smadl channel when the pathway is 
stimulated with TGF-P [25,27,28]. Human keratinocytes HaCaT 
cells exhibit a high basal concentration of Smad7 that is minimally 
affected upon treatment with TGF-P [29]. 

Robustness is a fundamental characteristic of biological systems, 
defining their ability to maintain normal function despite 
perturbations of their components [30]. In the TGF-P signal 
transduction pathway, negative feedback control in specific cases 
has been shown to confer robustness to the system, reducing 
phenotypic variability in cell populations [16], and to decrease the 
sensitivity of the signaling output to perturbations of its parameters 
[15]. A typical avenue to analyze the systemic robustness of a 
system is to perturb the model parameters, quantify this variation 
from the nominal parameter set, and assess the properties of the 
system output in the perturbed state. This type of approach has 
been used to study models of bacterial chemotaxis [31], the 
mitogen-activated protein kinase (MAPK) cascade [32], the 
interferon-gamma (IFN-y) induced Janus kinase (JAK) signal 
transducers and activators of transcription (STAT) pathway [33], 
the B-cell lymphoma 2 (Bcl-2) apoptotic switch [34], and the 
epidermal growth factor receptor (EGFR) signaling network [35]. 

Here we characterize all the potential Smad7-mediated negative 
feedback network motifs of the TGF-P signaling pathway and 
study their effects on the signaling dynamics, robustness, and 
sensitivity of a detailed mathematical model of the pathway. This 
type of study is notably important because network motifs, defined 
by a particular pattern of biochemical interactions, may reflect 
important functional properties of the system [36] . We investigate 
the dynamic response by exposing the system to two different 
extracellular ligand conditions, namely to stimulation with TGF-P 
ligand alone and to co-stimulation with TGF-P and BMP ligands. 
We do not consider Smad6 explicitly because its effects can be 
taken into account effectively by changing the strength of Smad7 
interactions and it does not give rise to any new negative feedback 
network motif. The robustness analysis considers a single measure 
for the whole system and provides insight into how the 
architecture of the network shapes the model's response to 
systemic parameter perturbation. Subsequently, to elucidate the 
effects of individual parameter perturbation on the model output, 
we use a global sensitivity analysis, which evaluates these effects 
within a large parameter space [37]. As the output of the model we 
focus on different properties of the nuclear concentration of 
phosphorylated R-Smad-Smad4 complexes, as they act as 
transcription factors to control the expression of the target genes. 
Our analysis provides a comprehensive examination into the 
effects of distinct negative feedback network motifs in defining the 
system behavior in the TGF-P signal transduction pathway. 

Methods 

TGF-p signal transduction pathway model 

We consider the detailed model developed in Ref. [25] to assess 
the Smad-dependent response to treatment with TGF-P and BMP 
ligands. We have shown elsewhere [25] that this detailed 
computational model accurately reproduces the diverse behavior 
of experimental datasets for human keratinocytes (HaCaT), bovine 



aortic endothelial cells (BAEC), and mouse mesenchymal cells 
(C2C12). The model includes three modules of the signaling 
pathway, namely receptor trafficking, nucleocytoplasmic shuttling 
of two parallel R-Smad channels, and a Smad7-based negative 
feedback loop. Signaling is initiated when a TGF-P or BMP ligand 
binds to its type II receptor, denoted as RII T or RII B , respectively. 
This ligand-receptor complex then recruits a type I receptor, 
either E1 1T , RI IB , or RI 2 . The former two type I receptors signal 
through the Smad 1 channel after binding to TGYfi-RIIf or BMP- 
Rlljj ligand-receptor complexes, respectively, while the latter 
signals through the Smad2 channel after binding to the TGFP- 
RII T complex. The resulting active heteromeric ligand-receptor 
complexes are denoted by Cjt, GYb, or C 2 , respectively, where the 
subscript is identical to that of the type I receptor within the 
complex. The active ligand-receptor complexes are then internal- 
ized into the early endosome, which provides a platform to 
efficiently phosphorylate cytosolic Smadl (Sl c ) or Smad2 (S2 C ). We 
use the subscripts c and n to indicate cytosolic and nuclear species, 
respectively, and the prefix p to denote phosphorylated species. In 
the cytosol, pSl c and pS2 c bind to Smad4 (S4 e ) to form the pSlS4 c 
and pS2S4 c complexes, respectively, which translocate into the 
nucleus. The complexes pSlS4 n and pS2S4 n then bind to DNA and 
activate the expression of Smad7 (ST), which irreversibly binds to 
surface ligand-receptor complexes (C ]T> C ]B , or C 2 ), preventing 
their association with and phosphorylation of R-Smad proteins, 
and targeting the active ligand-receptor complexes for degrada- 
tion. Therefore, the negative feedback loop is initiated with 
expression of Smad7 through the Smadl and Smad2 channels by 
way of pSlS4 n and pS2S4 m respectively. Smad7 proteins then 
inhibit the activation of the Smadl channel by binding to £y r and 
C 1B , while binding to C 2 inhibits Smad2 channel activation. A 
schematic representation of the model is shown in Figure 1 , where 
arrows indicate each modeled reaction in the pathway. Reactions 
are mathematically represented with mass-action kinetics, which 
are then used to form the system of ordinary differential equations 
(ODEs) to track the change in concentration of each modeled 
species over time (Table 1). 

To define the reference parameter set for the model, we use as 
starting point the parameter values defined for HaCaT cells in 
Ref. [25]. We assume that Smadl and Smad2 have the same 
abundance and reaction dynamics. For the association rate of 
Smad7 with C 2 (k 2 o a2 ), we use 1.50xl0 -4 molecules -1 min -1 , 
which corresponds to the initial estimated value before optimiza- 
tion for HaCaT cells [25]. This higher affinity value provides 
stronger inhibition through the negative feedback loop, but does 
not qualitatively affect the robustness and sensitivity results (see 
Text SI). We then set the parameters governing Smad7 expression 
by, and inhibition of, the Smadl channel, namely k 2 oa,iT> k 2 0a,iBi 
kjjpj, and K A1 to the corresponding values for the Smad2-channel- 
associated counterparts: k 20a ,2, hip,2> an d Ka,2- 

Prior to ligand stimulation, we determine the steady state 
solution for the system of ODEs (Table 1) by setting each time- 
derivative to zero and solving the linear system of equations that 
arises with TGFji and BMP equal to zero, the pre-stimulus 
conditions, using the 'linalg.solve' method in Numpy 1.6.2 (http:// 
numpy.scipy.org) with Python 2.7.3 (http://www.python.org). 
Upon adding the ligand, we numerically solve the system of 
ODEs using the CVODE method in the SUNDIALS 2.5.0 
package [38]. Thus, we focus on the typical experimental 
conditions that measure the response of the system to a sudden 
change of the ligand (TGF-P or TGF-P and BMP) concentration 
from zero to a saturating value (at time t = 0 hours) that is kept 
constant afterwards. 
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Figure 1. Schematic illustration of Smad-dependent TGF-p signal transduction pathway model. Arrows indicate reaction steps along the 
pathway and are labeled with the rate constant for the reaction. Dashed arrows denote Smad7 synthesis through gene expression. We use overbars 
to represent internalized receptor species in the endosome and the symbol "/" to indicate "or" in grouping the C 1T and C, 8 ligand-receptor 
complexes as C 1T/B . Different colors group the three modules, where blue indicates receptor trafficking, green indicates Smad nucleocytoplasmic 
shuttling, and red indicates negative feedback. 
doi:10.1371/journal.pone.0083531.g001 



Negative feedback network motifs 

We have considered nine distinct network motifs for the Smad7- 
based negative feedback loop in the TGF-P signaling pathway, 
which are schematically represented in Figure 2. These include 
three unbiased network motifs where inhibition equally affects 
both R-Smad channels, denoted here by "no degradation" (ND), 
"no feedback" (NF), and "total feedback" (TF). The ND network 
motif captures a system lacking Smad7, eliminating Smad7- 
dependent ligand-induced degradation of active ligand-receptor 
complexes from the model. In the NF network motif, Smad7 is 
kept at a constant concentration, analogous to a system with 
saturated levels of Smad7. Therefore, the NF network motif 



provides inhibition without the Smad7 negative feedback loop. For 
the TF network motif, Smad7 expression is activated by, and 
inhibits, both R-Smad channels, which corresponds to the 
complete network represented in Figure 1 . 

The remaining six network motifs include biased inhibition of a 
single R-Smad channel. Two of them, denoted by Al and A2, 
correspond to an auto-regulation negative feedback network motif, 
where Smad7 expression is activated by, and inhibits, the same R- 
Smad channel (Smadl for Al and Smad2 for A2). The C 1 and C2 
network motifs implement a cross-regulation negative feedback 
loop in which Smad7 expression is activated by the channel it does 
not inhibit. Finally, the network motifs denoted by SI and S2 
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Table 1. System of ODEs for each modeled species. 
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Equations track the change in concentration of each modeled species over time f. Overbars indicate receptor species internalized in the endosome. We prefix a species 
variable name with p to denote that it is phosphorylated and use the subscripts c and n represent cytoplasmic and nuclear species, respectively. 
doi:1 0.1 371/journal.pone.0083531 .tOOl 



correspond to a negative feedback loop with a single target, where 
Smad7 expression is activated by both R-Smad channels, but 
inhibits only one of the channels. In all the cases, RIjb and Rljr 
an not separated in the examined variations of feedback motifs 
because Smad7 is supposed to bind to both of them in the same 
way. Specifically, both of them have the same type of L45 loop on 
the kinase domain, which determines the interaction with Smad7 
through the L3 loop on the MH2 domain [39]. Since the kinetic 
parameters of Smad7 inhibition have not been fully characterized, 
based on this structural evidence, the most neutral assumption is to 
consider the same inhibitory kinetics for both of these type I 
receptors. 

In order to implement each specific network motif, we eliminate 
particular Smad7-based processes by setting their rate constant 
values to zero (Table 2). In addition, the initial concentration of 
Smad7 is set to about 2,633 molecules for the NF network motif. 
This leads to a degradation rate of 0.395 min 1 for the active 
ligand-receptor complexes, as in Ref. [13], for the reference values 
for k 20a jr, k 20a lB , and k 20a2 of Table 2. The initial Smad7 
concentration for the other eight network motifs is zero. We 
include a representation of the model built with CellDesigner 4.2 
[40] and the corresponding parameter values in File SI. This 
SBML file [41] corresponds to the TF network motif with its 
parameter set and initial conditions. The other cases are obtained 
by just setting the parameter values to those of Table 2. 

Robustness analysis 

The complexes pSlS4 n and pS2S4 n act as transcription factors to 
control the expression of hundreds of genes [5] . Thus, we focus on 
the different properties of the concentrations of these two species 
as the output of the model in our robustness analysis. Specifically, 
we investigate the following three properties: the peak species 



concentration (nip), the time at which the peak concentration is 
reached (ml), and the signal duration (mj) as defined in Ref. [42]. 
The peak species concentration, m p , is defined as 



m p = max(I(()), 



(1) 



which corresponds to the maximum of X(t), which is the 
concentration of pSlS4 n or pS2S4 n as a function of the simulation 
time t, with 0 < / < T. In the simulations, time t = 0 corresponds to 
the time at which the ligand is added and t = T corresponds to the 
total simulation time. We compute m, as 



--t(max(X(t))), 



(2) 



which corresponds to the time at which X(t) reaches its maximum 
value nip within the simulation time. The final metric, nij, is given 
by 



Jo r t 2 X(t)dt 



J 0 r tX(t)dt\ 2 
WW* 



(3) 



The signal duration metric captures the spread of X(t) about its 
peak [42] . In all three cases, we used a total simulation time T of 
24 hours. 

In order to assess the robustness of each network motif to 
parameter perturbation, we use the approach developed in Ref. 
[33] for the analysis of the IFN-y induced JAK-ST AT pathway. 
This method directly compares the parameter variation measure, 
defined by the logarithmic change in parameter values upon 
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Figure 2. Schematic representation of the nine negative feedback loop network motifs. Horizontal arrows represent mass flow into and 
out of a species. Vertical arrows and flat-head lines terminating on a horizontal arrow denote activation or repression, respectively, of the targeted 
process. Processes common to all network motifs are drawn in black, while the unique negative feedback loop processes are colored red. 
doi:10.1371/journal.pone.0083531.g002 



random perturbation, with the resulting metric output of the 
system. The parameter variation measure Af, is defined as 



7=1 



logic 

Pi 



(4) 



where n is the number of nonzero parameters, pjf is the /' 
parameter of the i ik test parameter set, and/y^ is the/* parameter 
in the reference parameter set [33]. For each parameter, we 
compute its test value by generating a random variable x with a 
uniform distribution within the range [— 1, 1] and multiplying the 
reference parameter value by 10*. This is an adaptation of the 
sampling method used in Ref. [33], where the authors generate * 
from a standard normal distribution. Our approach adapts the 
parameter space approach used in the parameter optimization 
routine from Ref. [25] . Specifically, we generate the test parameter 
set within a parameter space defined as ± 1 order of magnitude for 
each parameter, including the Smad7-associated parameters 



20a,lT, k'20a,lB, K 20a,2, K lip,h K lip,2: 



kupj, hip,2> K A h and K A J). We generate 



4,000 test parameter sets, as it is a sufficiently large number of 
samples for the results in this analysis to converge (see Text SI). 

For each test parameter set i, we simulate the model and define 
its output as C,. The different cases studied correspond to diverse 
ligand stimulation conditions, properties of the output, and output 
species for each of the nine negative feedback network motifs. This 
leads to a total of 108 cases resulting from all the possible 
combinations of taking one item from each of the following four 
groups for each of the 4,000 test parameter sets: 

(i) nine negative feedback network motifs: ND, NF, TF, Al, 
A2, CI, C2, SI, and S2; 

(ii) two ligand stimulation conditions: TGF-P alone and co- 
stimulation with TGF-fi and BMP together; 

(iii) two species: pSlS4 n and pS2S4„; 

(iv) three metrics: trip, m h and tttj. 

Therefore, the output C t represents one of three metrics for one 
of the two species for one of the two ligand stimulation conditions 
for one of the nine negative feedback network motifs, where the 
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Table 2. Parameter values of Smad7-related processes for the negative feedback network motifs. 







Parameter 


ND 


NF 


TF 


A1 


CI 


SI 


A2 


C2 


S2 


k20a.1T 

(molec mm ) 


0 


1.50x10 -4 


1.50x10 -4 


1.50X10 -4 


1.50 XlO -4 


1.50 xlO -4 


0 


0 


0 


^20a, IB 

(molec mm ) 


0 


1.50x10 -4 


1.50x10 -4 


1.50x10 -4 


1.50x10 -4 


1.50 xlO -4 


0 


0 


0 


l<20a,2 

(molec min ) 


0 


1.50x10 -4 


1.50x10 -4 


0 


0 


0 


1.50x10 -4 


1.50 xlO -4 


1.50 xlO -4 


hip, ; 

(molec min -1 ) 


0 


0 


8.53 xlO 3 


8.53 xlO 3 


0 


8.53 xlO 3 


0 


8.53 x10 s 


8.53 x10 s 


(molec min -1 ) 


0 


0 


8.53x10 3 


0 


8.53 x10 3 


8.53x10 3 


8.53 x10 s 


0 


8.53 x10 s 


(molec -1 ) 


0 


0 


1.03x10 -6 


1.03x10 -6 


0 


1.03 xlO -6 
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subscript i indicates one of the 4000 test parameter sets. One 
particular case consists of, for instance, the study of the ND 
network motif with just TGF-P stimulation and assessing as output 
the pSlS4 n species focusing on the m p metric for each test 
parameter set, yielding 4,000 values of C, for this case. 

For each case, we define a single robustness measure R as the 
ratio of the variance of all 4,000 values of A/, in Equation (4) with 
respect to that of Q. This is mathematically expressed as 



R 



E (Mj—M) 

/=i 

E(c,-c) 2 



(5) 



where JVis the number of parameter sets, equal to 4,000 here, and 
M and C are the mean values of A/, and C„ respectively. Larger 
values of R would indicate that the system is more robust to 
parameter perturbation, as the variance of the model output (C) is 
minimal compared to that of the parameter variation measure (A/). 
We have calculated this robustness measure for all of the 108 



Sensitivity analysis 

We have used a global sensitivity analysis to investigate how 
perturbations of individual parameters in the parameter space 
affect the model output [37]. Specifically, we performed a 
derivative-based global sensitivity analysis, which samples the 
effects of local parameter perturbation within the parameter space 
[43] . To estimate the effects of local parameter perturbation, we 
compute the scaled sensitivity coefficients [44] given by 



(6) 



where corresponds to the value of the/* parameter, k p of sample 
i with its corresponding model output C„ defined as in the 
robustness analysis. To approximate the partial derivative, we 
evaluate the model with 0.5 percent perturbations of each 
parameter h about its original value in sample i and calculate 
the finite central difference of the sensitivity metric [45] . 



For each parameter h, we compute the local scaled sensitivity 
coefficients E kj) for the 4,000 parameter values by randomly 
sampling its value in the parameter space, as described in the 
previous section for the robustness analysis, and use them to obtain 
three sensitivity measures [43]. The first measure, denoted by A k , 
averages the absolute value of the local sensitivity coefficients and 
is given by 



A k .= ~Y \E k ..\ 

K l TV t— 1 I J'' I 



(7) 



where jV is the number of parameter samples. The second 
measure, denoted by S k . , corresponds to the standard deviation of 

the local sensitivity coefficients and is given by 



*j ~ \ 



N ' 



£(W-<H- 



(3) 



The third measure, Gkj, corresponds to the sum of the squares 
of the first two measures and is given by 



(9) 



As the final measure combines the first two, we use G k , as the 
global sensitivity coefficient [43] . 

We implemented these robustness and sensitivity analyses using 
Python 2.7.3 (http://www.python.org), Numpy 1.6.2 (http:// 
numpy.scipy.org), and Scipy 0.10.0 [46]. 

Results 

Negative feedback network motifs' stimulation dynamics 

The specific implementation of the negative feedback loop 
differentially affects the dynamics of the model upon TGF-P 
stimulation (Figure 3). Network motifs with unbiased inhibition of 
the R-Smad channels, namely ND, NF, and TF, provide both 
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Figure 3. Model dynamics for stimulation with TGF-p. Simulation results upon single-ligand stimulation with TGF-p for the (A) ND, (B) NF, (C) 

TF, (D) A1, (E) CI, (F) SI, (G) A2, (H) C2, and (I) S2 network motifs. The simulation results of pS1S4 n , pS2S4 n , and S7 are displayed as black solid lines, 

black dashed lines, and red solid lines, respectively. 

doi:10.1371/journal.pone.0083531.g003 



Smadl and Smad2 channels with the same type of wiring upon 
single-ligand TGF-fi stimulation. Therefore, when the parameters 
for both channels are the same as in our case, the results for the 
dynamics of nuclear pSmadl-Smad4 and pSmad2-Smad4 species 
are identical (Figures 3A, 3B and 3C). In contrast, in the case of 
network motifs with biased inhibition, i.e., those where Smad7 
inhibits only one channel, such as Al, A2, CI, C2, SI, and S2, the 
system response varies between the two channels (Figures 3D, 3E, 
3F, 3G, 3H, and 31). In these network motifs, as Smad7 alternates 
its target between the Smadl and Smad2 channels, the dynamic 
response of one R-Smad channel mirrors that of the other. For 
example, in the case of the auto-regulatory negative feedback 
network motifs, the pSlS4„ response in the network motif A 1 is 
identical to that of pS2S4„ for the network motif A2 (Figures 3D 
and 3G). 

Upon co-stimulation of the system with both TGF-P and BMP 
ligands, the pSlS4 n and pS2S4„ signals are no longer identical for 
the ND, NF, and TF network motifs due to the additional 
activation of the Smadl channel by BMP (Figures 4A, 4B, and 
4C). In the case of the network motifs with biased architectures, 
namely Al, A2, CI, C2, SI, and S2, the system loses the mirrored 
behavior observed for R-Smad dynamics for stimulation with just 
the TGF-P ligand (Figures 4D, 4E, 4F, 4G, 4H, and 41). Under 
both stimulation conditions, just with TGF-P or with TGF-P and 
BMP, the signaling species targeted by the negative feedback loop 
exhibits a pronounced transient response, while the other species 
exhibits a more sustained response. 



Robustness analysis 

Our analysis reveals distinct variability in robustness among the 
negative feedback network motifs (Figure 5). When considering the 
peak species concentration (m f metric; Equation(l)) as the system 
output, the ND network motif consistently displays low values of 
the robustness measure R, indicating low robustness to parameter 
perturbation under the two stimulation conditions studied, namely 
stimulation with TGF-P only and co-stimulation with TGF-P and 
BMP ligands, for both pSlS4„ and pS2S4 n species (Figure 5A). The 
NF and TF network motifs lead to higher values of R, consistent 
with an increase of the system robustness compared to that of the 
ND one. Note that the R values are slightly different between the 
pSlS4n and the pS2S4n cases with the stimulation of just TGF-P in 
the case of NF and TF network motifs because they have been 
computed with different realizations of the random values of each 
of the parameters in the 4000 test parameter sets. As the number 
of test parameter sets goes to infinity, the difference between these 
values should vanish. 

The biased Al , A2, C 1 , C2, S 1 , and S2 network motifs exhibit a 
distinct pattern where the negative feedback loop increases the 
robustness of the species it targets in each case. For example, 
analyzing the pSlS4„ maximum concentration as the model output 
upon stimulation with TGF-P (Figure 5A, top row), our results 
indicate that Al, CI, and SI network motifs are more robust than 
those where the negative feedback loop inhibits the Smad2 
channel (A2, C2, and S2) and vice versa. Stimulating the system with 
both TGF-P and BMP ligands results in lower values of the 
robustness measure for the Smadl channel (pSlS4„ species) for all 
network motifs when compared to stimulation with TGF-P alone, 
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Figure 4. Model dynamics for stimulation with TGF-p and BMP. Simulation results upon co-stimulation with TGF-p and BMP for the (A) ND, 

(B) NF, (C) TF, (D) A1, (E) CI, (F) SI, (G) A2, (H) C2, and (I) S2 network motifs. The simulation results of pS1S4 n , pS2S4 n , and S7 are displayed as black 

solid lines, black dashed lines, and red solid lines, respectively. 
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while the robustness of the Smad2 channel (fiS2S4„ species) 
increases, but is less affected by co-stimulation. 

Considering other metrics, specifically the time at which the 
maximum concentration for the different species is reached (m t 
metric; Equation(2)), results in a qualitatively similar system 
robustness when compared to the nip metric for all network motifs 
and stimulation conditions (Figure 5B). In the case of the signal 
duration (m d metric; Equation(3)), the robustness analysis reveals a 
decreased robustness for the negative feedback loop, in contrast to 
the behavior observed for the nip and m, metrics (Figure 5C). In this 
case, the NF and ND network motifs display the highest robustness 
values, while robustness decreases for the network motifs in which 
the negative feedback loop inhibits the evaluated channel. 

Sensitivity analysis 

The results of the sensitivity analysis are shown in Figure 6, 
Figure SI, and Figure S2 for the nip, m t , and m d metrics, 
respectively. For the peak species concentration metric nip, our 
analysis reveals a high network-motif-independent sensitivity, as 
given by a high value of the global sensitivity coefficients (Gfy 
measure; Equation (9)), for multiple parameters, including the 
synthesis and degradation rate constants of the R-Smads (k syni Rs 
and k desRS ) and Smad4 {k m> s4 an d k degi s4) (Figure 6). Additionally, 
the model is sensitive to the parameters governing the reversible 
association of nuclear phosphorylated R-Smad with Smad4 (kio a 
and kj 0 J), nucleocytoplasmic shutding of Smad4 (k ]4lm p and k 14ex p), 
and degradation of nuclear phosphorylated R-Smads (k J5deg ) for all 
network motifs. The model's sensitivity for the Smad7-related 
parameter values, namely k 20aJT , k 2llaJB , k 20<2 , k, p j, k ap:2 , K AJ , 



and K A2 , is motif-dependent. Specifically, the model is more 
sensitive to these processes if the negative feedback loop inhibits 
the evaluated channel. For example, when stimulating the system 
with TGF-P alone and assessing the pSlS4 n species as the output, 
this set of parameters has a higher sensitivity coefficient in the 
network motifs Al , CI, and S 1 than in the network motifs A2, C2, 
and S2. 

For the metric m, that characterizes the time at which the peak 
concentration is reached, the sensitivity analysis shows minimal 
differences in the sensitivity coefficients among all parameters 
(Figure SI). Nevertheless, it reveals a higher sensitivity for several 
groups of parameters in the network motifs where Smad7 does not 
inhibit activation of the evaluated species (e.g. A2, C2, and S2 
network motifs are more sensitive than Al, CI, and SI network 
motifs to perturbations of k 10a , k I0d , and k 15deg parameters when 
evaluating the pSlS4„ species). This effect is consistent with the 
results of the robustness analysis, in the sense that increases in the 
robustness measure R correspond with a decreased sensitivity 
coefficient. 

The sensitivity analysis for the signal duration metric m d 
indicates that the model exhibits a higher sensitivity to parameter 
perturbation when the negative feedback loop inhibits the 
evaluated species (Figure S2), which is consistent with the 
robustness analysis results for this metric. This effect is most 
significant for the Smad7-related parameters k 20aIT , k 20a lB , k 20a 2 , 
ki,pj, k lip 2 , K A i, and K A 2 . 
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Figure 5. Robustness analysis. Each negative feedback network 
motif is assessed for its robustness by computing the value of the 
robustness measure ft (Equation(5)) upon stimulation with TGF-p alone 
or together with BMP. To represent the model output, we use the (A) 
peak species concentration (m p metric; Equation(1)), (B) time of the 
peak species concentration (m t metric; Equation(2)), and (C) signal 
duration (m d metric; Equation(3)) as metrics of the pS1S4 n and pS2S4 n 
signal response. 
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Discussion 

The TGF-fi signal transduction pathway is extensively regulated 
to effectively control multiple cellular responses [2,6]. Here, we 
have characterized the effects of a key regulatory mechanism that 
controls the potential attenuation of the signal through a Smad7- 
dependent negative feedback loop. Our results show that 
variations in the network design of this negative feedback loop 
result in distinct dynamic behavior and response to parameter 
perturbation. Specifically, repression by the negative feedback 
loop results in a transient signal response, while the absence of this 
feedback results in a more sustained behavior. As regulation of 
gene expression by pSlS4 n and pS2S4 n is linked to this signaling 
property, requiring a prolonged signal to maximally activate 
transcription [47], negative feedback may be a central control 



mechanism for determining the long-term cellular response to 
ligand stimulation. 

The robustness analysis for the mp and m t metrics revealed that 
the presence of the negative feedback loop correlated with 
increased robustness for the signaling channel it represses. 
However, our results also show that the three signaling metrics 
we examined do not respond identically to parameter perturba- 
tion. In particular, the robustness of the signal duration metric 
decreases for the species repressed by the negative feedback loop. 
With negative feedback, the repressed signal exhibits a transient 
response with shorter signal duration than the sustained response 
of the non-repressed signal. This type of transient response is 
characteristic of cancer cells [48]. Our results suggest that the 
shorter signal duration is more significandy affected by parameter 
perturbation than the greater signal duration exhibited by the 
species without negative feedback repression. Together, the three 
metrics show that the robustness of the system to parameter 
perturbation is dependent upon the specific negative feedback 
network motif, although in a different manner for the mp and m l 
metrics than the m d metric. 

Similar to the robustness analysis, the sensitivity analysis shows 
varied results for the three signal metrics. Assessing the peak 
species concentration (nip metric) shows an increased sensitivity 
coefficient for the Smad7-based processes when the negative 
feedback loop represses the evaluated species. However, the 
robustness analysis for this metric shows an increased robustness 
with the addition of a negative feedback loop, revealing that, while 
the negative feedback loop itself may be sensitive to perturbation, 
its impact on the complete system results in additional robustness. 
With the m t metric, the specific negative feedback network motif 
has a minimal effect on the sensitivity coefficients for the majority 
of parameters. However, with several parameters (k^ nR& kj egiJis , 
k 10a , k ]0lt , and k 15deg ), the sensitivity analysis shows as well that the 
negative feedback loop establishes resistance to perturbations by 
lowering the sensitivity coefficient compared with the network 
motifs lacking repression of the evaluated species. 

These results, as a whole, show that the robustness and 
sensitivity properties depend both on the specific network motif 
and on the signal property of interest. This type of behavior is 
typically observed in other systems. For instance, in the classical 
example of bacterial chemotaxis [49] , it was observed that steady- 
state behavior and adaptation time are not robust, while the 
precision of adaptation is robust to changes in protein concentra- 
tions. In our case, as in the case of bacterial chemotaxis, the 
adaptation time, or signal duration, is not a robust property. The 
general rule is that an increase in robustness against some 
perturbations will be counterbalanced by a decrease of robustness 
against other perturbations [50]. This rule is epitomized by the 
Bode integral formula, which represents conservation of sensitivity 
of a negative feedback system along the frequency axis [51,52]. 

The available experimental data indicates that the negative 
feedback loop exists in several forms in different cell lines. In 
particular, bovine aortic endothelial cells exhibit an auto- 
regulatory negative feedback loop in which Smad7 is expressed 
through activation of the Smadl channel and inhibits further 
activation of the same R-Smad channel [26]. In the mouse 
myoblast ceU line C2C12, TGF-P does not activate BMP- 
responsive reporter genes [28], suggesting that the ligand induces 
transcriptional activity through the Smad2 channel. Results from 
experiments tracking the dynamics of R-Smad phosphorylation 
[27] and from computational modeling of the signaling behavior 
in these cells [25] suggest that Smad7 primarily inhibits activation 
of the Smadl channel. Together, these findings provide evidence 
that C2C12 cells possess a cross-regulatory negative feedback loop 
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Figure 6. Sensitivity analysis for the peak species response metric nip. We compute the sensitivity coefficient Gj ( (Equation (9)) for each 
model parameter for the nine negative feedback network motifs upon stimulation with TGF-p alone or together with BMP. We use the metric m p 
(Equation(l)) to assess the signal response of pS1S4 n and pS2S4„. Empty spaces (white rectangles) indicate a sensitivity coefficient equal to zero or a 
parameter value set to zero from Table 2. 
doi:10.1371/journal.pone.0083531.g006 



in which Smad7 is expressed through the Smad2 channel, but 
inhibits the Smadl channel when the pathway is activated with 
TGF-p. In contrast, human keratinocyte HaCaT cells exhibit a 
high basal concentration of Smad7 that is minimally affected upon 
treatment with TGF-P [29]. 

Our results are important because many components of the 
TGF-P signaling pathway are mutated, downregulated, or 
overexpressed in multiple diseases, such as the TGF-P receptors, 
R-Smads, Smad4, and Smad7 proteins in a variety of cancer types 
[3]. The sensitivity analysis captures the effects of these 
perturbations, quantifying how the model responds to variations 
in the pathway reactions. Indeed, our results from the sensitivity 
analysis identify several processes with high sensitivity coefficients, 
which are often dysregulated in cancer cells. For instance, missense 
mutations in the Smad4 gene found in pancreatic cancer cells are 



associated with reduced nuclear translocation [53]. The model 
describes nuclear translocation of Smad4 with the rate constants 
kj4, m f an d kj4 a p, both of which display among the highest 
sensitivity coefficients with the nip and m d metrics. Furthermore, 
missense mutations in the Smad2 and Smad4 genes occurring in 
colon and pancreatic cancer cells, respectively, have been reported 
to inhibit association of Smad2 with Smad4 [54] . The sensitivity 
analysis results show high sensitivity coefficients for the rate 
constants governing this process as well, where k 10a and k 10d 
regulate the nuclear association of phosphorylated Smad2 and 
Smad4. Notably, the sensitivity coefficients for k 10a and k 10d arc 
dependent on the specific negative feedback network motif, which 
is most significantly observed with the m t metric, where network 
motifs lacking inhibition of the evaluated species display a higher 
sensitivity coefficient than those where the negative feedback loop 
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inhibits the evaluated species. This correlation between model 
sensitivity and pathway mutation indicates that our analysis may 
be used to elucidate which processes are involved in the transition 
from normal to pathological states in a variety of cell types that 
exhibit the different negative feedback network motifs. 

The results of the sensitivity analysis additionally provide a tool 
for determining novel targets in the pathway for therapeutic 
intervention. Potential therapeutic targets are defined as those 
where perturbations significandy affect the signaling response, 
such that administering treatment will maximally impact the 
dynamic behavior. By applying this analysis to the different 
negative feedback network motifs our results can be used to 
identify the therapeutic potential for targeting processes in a 
variety of cell types. 

Supporting Information 

Figure SI 

Sensitivity analysis for the time of the peak species 
response metric m t . We perform the same analysis described 
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